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Abstract 

Measurements of the velocity field created by a shallow bump on a wall revealed that an 
energy peak in the spanwise spectrum associated with the driver decays and an initially small- 
amplitude secondary mode rapidly grows with distance downstream of the bump. Linear 
theories could not provide an explanation for this growing mode. The present Navier-Stokes 
simulation replicates and confirms the experimental results. Insight into the structure of the 
flow was obtained from a study of the results of the calculations and is presented. 
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I. Introduction 


In order to determine whether an eigenmode of the linearized boundary-layer equations 
or an Orr-Sommerfeld approach was appropriate for studying the disturbance field created 
by a localized boundary perturbation of boundary-layer flow, Gaster, Grosch, and Jackson 1 
performed an experiment to find out what the solution should look like. Surprisingly, the 
experimental results did not seem to be consistent with either model. A significant result of 
the experiment was that a primary spanwise mode associated with the diaphragm dimension 
decayed slightly with downstream distance and a small secondary mode rapidly grew with 
downstream distance. To rule out any possible anomaly in the experiment, the present 
computational study was initiated to duplicate the experimental results. Equally important, 
the results of the computation permit a detailed study of the flow field structure and yield 
insight into the physics of the flow. The computations involve solving the unsteady nonlinear 
Navier-Stokes equations for a spatially-growing boundary-layer flow and should be equivalent 
to an ideal experimental study. 


II. Overview of Experimental Conditions 


Although a detailed description of the experimental conditions was given by Gaster et 
al. 1 , a synopsis of the experimental parameters important for the computations is provided 
here. In the experiment, the bump was located 400 mm from the leading edge of the 
flat plate. At this location, the boundary and displacement thicknesses of the undisturbed 
flow near the bump were 8 — 2.88 mm and 8* — 0.99 mm, respectively. The associated 
boundary-layer Reynolds numbers near the bump were Rs = 3480 and Rs- = 1196. A 
silicon rubber diaphragm of 20 mm (~ 20 J*) was used to force the disturbance. Most of 
the diaphragm motion occurred over 10-15 mm of the center. The amplitude of the bump 
motion was about 0.1 mm (100 /mi), which is a typical height of the roughness element used 
in receptivity experiments (see Saric 2 , section 3.1.1). A measure of the disturbance amplitude 
near the bump was predicted to be about u = 4.9% of the freestream velocity. Although a 
stationary bump would be preferable, a forcing frequency of 2 Hz was used to discriminate 
the signal created by the bump from the background noise present to some degree in all 
experiments. There are, of course, Tollmien-Schlichting modes with a 2 Hz frequency but 
these are highly damped at the Reynolds numbers of the experiment. The 2 Hz frequency is 
well below that of any growing Tollmien-Schlichting mode. The 4.9% disturbance amplitude 
is too large to enable a comparison of the experimental disturbance with receptivity theory 
which, to date, is based on the infinitesimal small-amplitude assumption. It is possible that 
a sufficiently large disturbance at very low frequency or even a steady disturbance could 
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cause a by-pass transition, but no evidence of transition or turbulence was observed in the 
experiment or in the computations reported here. Measurement stations were setup about 70 
and 105 boundary-layer thicknesses downstream of the bump, (~ 200 and 300 displacement 
thicknesses). These are station B, section b-b and station A, section c-c, respectively, as 
shown in figure 1 of Gaster et al 1 . The Reynolds numbers at these measuring stations 
were approximately R$* = 1466 and R$* — 1585. Detailed measurements of the streamwise 
component profiles in both the normal and spanwise directions were made at these stations. 
Less detailed measurements were made over a larger area. 


III. Numerical Method of Solution 


The numerical techniques required for the simulation and the disturbance forcing are 
briefly discussed in this section. For a detailed description of the spatial DNS (Navier- Stokes) 
approach used for this study, refer to Joslin, Streett and Chang. 3,4 The instantaneous ve- 
locities u — ( u , v , w ) and the pressure p are decomposed into steady base and disturbance 
components. The base flow is given by velocities U_ = (C, V, W) and the pressure P‘, the dis- 
turbance is given by velocities u = (u, v, w ) and the pressure p. The velocities correspond to 
the coordinate system x — ( x , y , z), where x is the streamwise direction, y is the wall-normal 
direction, and z is the spanwise direction. The base flow for the flat plate can be reason- 
ably approximated by the Blasius similarity solution U_ = ( U , V, 0) and the disturbance flow 
is found by solving the three-dimensional, incompressible Navier-Stokes equations. These 
equations are the momentum equations 

u t + {u- V)u + (U • V)u + ( M • V)C = -Vp + i V 2 n (1) 

and the continuity equation 

V • u = 0 (2) 

The boundary conditions in the far-field are 

u — > 0 as y -> oo (3a) 


and the conditions at the wall are 

u — u c at y — 0 (3 b) 

where u c = 0, except for the portion of the wall which models the bump. The Reynolds 
number R — is based on the boundary-layer displacement thickness at the inflow of 

the computational domain, the freestream velocity Coo and the kinematic viscosity v. 

To solve equations (1-3) computationally, the spatial discretization entails a Chebyshev 
collocation grid in the wall-normal direction, fourth-order finite differences for the pressure 
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equation, sixth-order compact differences for the momentum equations in the streamwise 
direction, and a Fourier sine and cosine series in the spanwise direction on a staggered 
grid. 4 For time marching, a time-splitting procedure is used with implicit Crank-Nicolson 
differencing for normal diffusion terms and an explicit three-stage Runge-Kutta method. 5 The 
influence-matrix technique is employed to solve the resulting pressure equation (Helmholtz- 
Neumann problem). 6,7 At the inflow boundary, the mean base flow is forced and, at the 
outflow, the buffer-domain technique of Streett and Macaraeg 8 is used. 

The boundary conditions that have to be imposed to represent an oscillating bump are 
a streamwise velocity perturbation related to the mean shear and the bump height together 
with the normal velocity of the bump as described in Gaster et al. 1 Although the experiments 
had to use a low-frequency oscillating bump, the computations can use a stationary bump. 
Hence, the boundary conditions reduce to 

w(a:,0, z) — — h(x,z)— (4) 

The form h(x, z ) = v w sin(«T) 3 sin(^) 3 is imposed for the bump shape, which yields a compu- 
tationally smooth bump. The bump height is given by the amplitude v w = 10%. Although 
no attempt was made to exactly match the ingested disturbance amplitude of 4.9% in the 
experiment, an amplitude of approximately 3.4% was observed in the computations. The 
streamwise length of the bump was 15.95* and the spanwise half-length was 6.55*. As will 
be seen in the results section of this paper, the somewhat arbitrary selection of the ampli- 
tude and shape of the computational bump did not have an adverse effect on the desired 
comparison with experiments. 

Figure 1 is a sketch of the computational domain. In the experiment the displacement 
thickness at the bump, 5* was approximately 1 mm and the boundary layer Reynolds number 
at that location, R = was 1200. We choose 5* as the length scale for the computation 

and set the Reynolds number in equation (1) to be 1200. As shown in figure 1, the far-field 
boundary was located 505* from the wall, the streamwise extent of the domain was 5005* 
from the inflow, and the spanwise extent of the domain was 255*. This spanwise extent is 
shown by the dotted line in figure 1. Along this surface a symmetry condition on the flow 
field was applied. Thus the effective spanwise extent of the computational domain was 505*, 
as shown in figure 1. The center of the bump was positioned on the symmetry boundary 
at 40.95* from the inflow, again as shown in figure 1. With the computational scale chosen 
to be 5* = 1 mm, all of the dimensions shown in figure 1 of Gaster et al. 1 can be directly 
translated into 5* units. The lines labeled a - a, b - b and c - c in figure 1 are the similarly 
labeled lines shown in figure 1 of Gaster et al. 1 . 

The choice of grid, computational domain size, and time-step size were based on pre- 
vious experience described in Joslin, Streett and Chang 3,4 for unsteady disturbances and 
in Joslin and Streett 9 for a stationary disturbance. The simulation used a coarse grid of 
661 streamwise, 61 wall-normal, and 20 spanwise grid points (spanwise symmetric). For the 
time marching, a time-step size of 0.2 is chosen for the three-stage Runge-Kutta method. 
The coarse grid computation required 44 Cray 2 hrs with a single processor to converge to 


3 



a time-independent solution. In addition to the coarse grid calculation, a grid refinement 
simulation was performed to verify the quantitative accuracy of the results of the course 
grid computations reported below. This second simulation was conducted with a grid of 
1321 streamwise, 81 wall-normal and 39 spanwise grid points; This translates into doubling 
the grid in the streamwise and spanwise directions; in the wall-normal direction Chebyshev 
series are used, which have coefficients that converge exponentially. 

Because the disturbance excitation is steady and the resulting disturbance modes are 
stationary, the fine-grid simulation had initial conditions which correspond to the course grid 
final results. If the course and fine grid results were time-independent and quantitatively 
similar, then significant computational savings (approximately 300 Cray Y/MP hours) can 
be realized with this choice of initial conditions for the fine-grid simulation. The fine-grid 
simulation was marched in time and the results were compared after 180 and 420 time steps. 
The results were identical, indicating that the fine grid simulation had converged after only 
180 time steps. The cost of this simulation was 19 Cray Y /MP hrs. As will be shown below 
in the Results section, the results of the coarse and fine grid computations were essentially 
identical. 


IV. Results 


The results shown in figure 2 are spanwise profiles of the streamwise velocity component 
in the spanwise direction at y ~ <£*, which is approximately the distance from the wall 
used in the experiments. This top view would have the bump placed at the bottom of the 
figure, and the flow direction is from bottom to top. As expected, there is a velocity deficit 
directly downstream of the bump and lobes of enhanced velocity on both sides of the bump. 
Note that the intensity of this deficit and lobes is decreasing with distance downstream of 
the bump. This qualitative picture matches the experimental observation, except there was 
some asymmetry in the experiments. 

Figure 3a shows the variation with downstream distance of the total energy of the dis- 
turbance generated by the bump, as obtained from the fine and course grid simulations. 
Clearly, quantitative agreement is observed (note that the ordinate has a logarithmic scale). 
Figure 3b shows the variation of the total energy and the square of the velocity components 
with downstream distance. The total energy is decreasing with distance downstream and 
the streamwise velocity component is clearly dominant compared with the insignificant wall- 
normal and spanwise components. In the experiments, only the streamwise velocity compo- 
nent was recorded and the discussion and conclusions of the flow were described based on 
the streamwise velocity. The computations clearly show that it is unnecessary to consider 
the wall-normal and spanwise velocities. 

Figure 4 shows the low-wavenumber modal decomposition of the streamwise velocity 
component in the spanwise direction. Confirming the experiments, the low-wavenumber 
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modes are growing with downstream distance; all other high- wavenumber modes (> 4) are 
decaying everywhere. The low-wavenumber modes are marked by a region of rapid growth 
followed by either an asymptote or decay region beyond the present computational domain. 
A growth in the (3 = 2 mode by a factor of 8 in magnitude was noted by Gaster et al. 1 Here, 
the dominant (3 = 2 mode has grown by over a factor of 6 in magnitude and has not reached 
its maximum value within the computational domain. 

The (3 = 2 velocity profiles, obtained from both the coarse and fine grid simulations, 
are shown with distance from the wall in figure 5 at various downstream distances. The 
results of both simulations are in excellent quantitative agreement and both simulations 
show modal growth consistent with the experiments. The profiles at R = 1576 and 1617 
qualitatively match the experiments in shape and have their peak near y ~ 25*, as do the 
experimental results (see figure 3 in Gaster et al 1 ). The magnitudes are, however, different. 
The results of the calculations shown in figure 5 have a peak value of about 5 x 10“ 5 , while 
the measurements show a peak value of 2 x 10~ 3 . 

The three dimensional structure of the flow field can be inferred from the results presented 
in figures 6, 7 and 8. All of these results are on the y-z plane at x = 503 which is slightly 
downstream of the section a-a as shown in figure 1 of Gaster et al 1 and of this paper. At this 
location R s * ~ 1388. The results of the calculation are obtained on a Chebyshev collocation 
grid in the wall-normal (y) direction. In the spanwise direction the computational results 
from z = 0 to 2 = 25 5* (with a symmetry boundary condition at z = 0) were “folded” about 
z = 0 in order to obtain the flow field in -25 5* <z< 255*. Although the far-field boundary 
is located at y = 505*, results are shown in 0 < y < 55* because the disturbance field is 
essentially confined to the boundary layer. Because it is more convenient in presenting the 
data, the computational results were interpolated onto a uniform grid in the y direction. It 
should be noted that figures 6, 7 and 8 are distorted by an, approximately, 10 to 1 stretching 
in the y direction as compared to the z direction. 

Figure 6 contains contours of u, the streamwise component of the velocity. In this plane, 
-1.69 x 10 -3 < u < 0.48 x 10 -3 and the minimum occurs on the centerline at y = 0.73, with 
the maxima being located at y = 0.73 and 2 = ±3.95. This is the same structure seen in the 
profiles of u at y = 5* shown in figure 2. As seen from figure 6, the streamwise component 
of the disturbance field is an upstream flowing “jet” on the centerline with downstream 
counter-flowing jets on both sides. The entire field is essentially confined to the inner part 
of the boundary layer (y < 2). 

The (v,w) vectors in the same x plane are shown in figure 7. The maximum of \/v 2 + w 2 
is 2.62 x 10~ 5 and occurs at y = 0.35 and z = ± 5.26. There is an inflow towards the region 
of the upstream “jet” along the centerline and an outflow between the downstream “jets” 
and the wall. Both the inflow and outflow are rather small compared to u; less than 1% 
of the maximum of u but extend over a very large region in the y-z plane. There is even 
a small, but appreciable, inflow at the top of the boundary layer. It should be noted that 
the cross stream flow is rather weak and the crossflow Reynolds number (see Saric 2 ) is very 
small. As can be seen from this figure, the boundary layer thickness of the cross flow (5c 
is of the same size as the boundary layer thickness (5) of the streamwise flow. In contrast 
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the magnitude of the crossflow velocity component (U c ) is very small compared to that of 
the mean flow (U 0 ). It is clear that the cross flow Reynolds number can be calculated by 
R c = ( 6 C /S)(U C /U 0 )R ■ This gives a value 0 ( 0.1 ). 

Contours of the streamwise (x) component of the vorticity, u> x , are plotted in figure 8. 
These were obtained by numerically differentiating v and w using a second order scheme. 
No smoothing of the results was done. It might have been expected that the numerical 
differentiation would induce substantial “noise”, but none is apparent in the results shown 
in figure 8. The maximum and minimum of u x are ± 1.64 x 10 -4 and lie at y = 0.70 and 
z — =p 3.95. Positive u> x indicates clockwise rotation and negative counter-clockwise. It is 
seen from the structures shown in this figure that the bump generates a pair of counter- 
rotating vortices just above and on either side of it. These pump fluid down towards the 
wall and into the upstream flowing “jet” of the disturbance field. Just above the main pair 
of vortices and slightly towards the centerline there are a weaker pair of oppositely rotating 
vortices. Between the main pair of vortices and the wall there is region of high vorticity due 
to the relatively strong outflow in the ± z directions. 

This basic structure of the flow field persists further downstream but is considerably 
weaker. This can be seen from the results shown in figures 9, 10 and 11. These results are on 
the y-z plane at x = 709 which is slightly downstream of the section c-c as shown in figure 
1 of Gaster et al 1 and of this paper. At this location Rs * « 1585. The contours of u are 
shown in figure 9 with the same contour levels as in figure 6. The disturbance u has the same 
general structure as at x = 503. However it is considerably weakened with the minimum 
of the upstream flowing “jet” only -0.63 x 10“ 3 . The maxima of the downstream jets are, 
however, slightly larger than at x = 503. The y position of the center of these “jets” is 1.06 
at this location as compared to 0.73 at x = 503 and the “jets” have diffused in the wall 
normal direction. The (v,w) vectors plotted in figure 10 also show uplift and spreading in y 
as well as a general weakening. The maximum of \/v 2 + w 2 at this x location is 0.38 x 10 — 5 , 
nearly seven times smaller than at x = 503. The z location of these maxima is exactly the 
same as at x = 503 but the y location is 0.62 as compared to y — 0.35 ar x = 503. Finally, 
figure 11 shows contours of u x at x = 709. The contour levels in this figure are one-tenth 
of those in figure 8. The general decrease and diffusive spreading in the vorticity is readily 
apparent. The main vortices have lifted further from the wall. Their centers are at the same 
z position as at x = 503 but the y position is now 1.20 instead of 0.70 and are at the same 
height as the secondary pair. The vorticity at the wall is also considerably weakened. 


V. Concluding Remarks 


The spatial evolution of the disturbance velocity field initiated from a shallow bump on 
a wall in a laminar boundary layer was computed by direct numerical simulation of the 


6 



incompressible Navier-Stokes equations. Comparison of results from coarse and fine grid 
simulations showed that the coarse grid simulation had converged. The evolution pattern 
and modal growth and decay trends were shown to be consistent with the experimental 
results of Caster et al. 1 

The three dimensional structure of the flow field was inferred by examining the velocity 
component and the streamwise component of the vorticity on two y-z planes downstream 
of the bump. It was seen that the bump generates, at least in the far field, a pair of 
counter-rotating vortices just above the wall and on either side of the bump location. These 
pump fluid down towards the wall and into an upstream flowing "jet ” of the disturbance 
field. Outside of this main “jet” there are a pair of weaker downstream flowing “jets”. Just 
above the main pair of vortices and slightly towards the centerline there are a weaker pair of 
oppositely rotating vortices. Between the main pair of vortices and the wall there is region 
of vorticity due to the outflow in the ± z directions. As this flow evolves downstream, the 
vortices lift from the wall, diffuse and weaken while maintaining their basic structure. The 
“jets” also weaken, diffuse and lift, as must happen because they are, in a sense, both cause 
and effect of the vortices. 

A theoretical study should be conducted to complete the understanding of this proposed 
linear transfer of energy between various spanwise modes. In addition, wind-tunnel and 
computational experiments should be conducted in order to understand the interaction of 
Tollmien-Schlichting waves with the bump-induced vorticity field. 
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Figure 1: Sketch of the computational domain showing its size and the location of the bump. 
The lines a- a, b-b and c - c show the location of the similarly labeled lines of Gaster et 
al. 1 , along which measurements were made. We present results of the simulation on planes 
including lines a — a and c — c. 




Figure 2: Streamwise velocity profiles in the streamwise/spanwise plane downstream of the 
bump. The profiles are at a height of y = S*. 
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Figure 3: Variation of the total disturbance energy with downstream distance as obtained 
from the coarse and fine grid computations. 
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Figure 4: The streamwise velocity component at three different downstream locations. The 
solid curves are the result of the coarse grid simulation and the dashed curves are the result 
of the fine grid simulation. These curves coincide almost everywhere in this figure. 
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Figure 5: Contours of u, the streamwise component of the disturbance velocity, on the y-z 
plane at x = 503. Contours with positive values of u are solid and those with negative values 
are dashed. The contours values are -1.6 x 10 -3 to -0.2 x 10 -3 in steps of 0.2 x 10“ 3 and 
from -0.1 x 10 -3 to 0.5 x 10 -3 in steps of 0.1 x 10~ 3 , excluding 0.0. The minimum value of 
u, -1.69 x 10~ 3 , occurs on the centerline at y = 0.73, with the maxima, 0.48 x 10 -3 , being 
located at y = 0.73 and z = ± 3.95. 
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Figure 6: Contours of u, the streamwise component of the disturbance velocity, on the y-z 
plane at x = 503. Contours with positive values of u are solid and those with negative values 
are dashed. The contours values are -1.6 x 10“ 3 to -0.2 x 10~ 3 in steps of 0.2 x 10 -3 and 
from -0.1 x 10~ 3 to 0.5 x 10~ 3 in steps of 0.1 x 10 -3 , excluding 0.0. The minimum value of 
u, -1.69 x 10 -3 , occurs on the centerline at y = 0.73, with the maxima, 0.48 x 10 -3 , being 
located at y = 0.73 and z = ± 3.95. 
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Figure 7: Vectors of (v,w) on the y-z plane at x = 503. The maxima of \/u 2 + tc 2 are located 
at y = 0.35 and z = ± 5.26 and have the value 2.62 xlO -5 . 
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Figure 8: Contours of u> x , the streamwise component of the disturbance vorticity, on the 
y-z plane at x = 503. Contours with positive values of u x are solid and those with negative 
values are dashed. Positive to x indicates clockwise rotation and negative counter-clockwise. 
The contours values are -1.6 x 10" 4 to 1.6 x 10~ 4 in steps of 0.2 x 10~ 4 , excluding 0.0. 
The minimum value of w*, -1.64 x 10" 4 , is at y = 0.70, and z = 3.95 and the maximum, 
1.64 x 10 -4 , is at y = 0.70 and 2 = -3.95. 
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Figure 9: Contours of u, the streamwise component of the disturbance velocity, on the y-z 
plane at x = 709. Contours with positive values of u are solid and those with negative values 
are dashed. The contours values are —1.6 x 10 -3 to —0.2 x 10~ 3 in steps of 0.2 x 10 -3 and 
from -0.1 x 10~ 3 to 0.5 x 10“ 3 in steps of 0.1 x 10 -3 , excluding 0.0. The minimum value of 
u, -0.63 x 10~ 3 , occurs on the centerline at y = 1.06, with the maxima, 0.51 x 10~ 3 , being 
located at y = 1.06 and z = ± 2.63. 
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Figure 10: Vectors of (v,w) on the y-z plane at x = 709. The maxima of a/u 2 + w 2 are 
located at y = 0.62 and z = ± 5.26 and have the value 0.38 x 10“ 5 . 





Figure 11: Contours of u> x , the streamwise component of the disturbance vorticity, on the y-z 
plane at x = 709. Contours with positive values of u are solid and those with negative values 
are dashed. Positive <x> x indicates clockwise rotation and negative counter-clockwise. The 
contour values are -1.6 x 10 -5 to 1.6 x 10~ 5 in steps of 0.2 x 10 -5 , excluding 0.0. Note that 
these contour levels are one tenth those of the contour plot of u x at x = 503. The minimum 
value of u> x , -1.87 x 10 -5 , is at y = 1.20, and z = 3.95 and the maximum, 1.87 x 10~ 5 , is at 
y = 1.20 and z = —3.95. 20 
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